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ABSTRACT 


The design of control laws to damp flexible structural modes requires accurate math 
models. Unlike the design of control laws for rigid body motion (e.g., where robust 
control is used to compensate for modeling inaccuracies), structural mode damping 
usually employs narrow band notch Alters. In order to obtain the required accuracy in 
the math model, maximum likelihood estimation technique Is employed to Improve 
the accuracy of the math model using flight data. This paper presents all phases of this 
methodology: (1) pre-flight analysis (i.e.. optimal Input signal design for flight test, 
sensor location determination, model reduction technique, etc.). (2) data collection and 
preprocessing, and (3) post-flight analysis (l.e., estimation technique and model 
verification). In addition, a discussion is presented herein of the software tools used for 
this study and the need for future study In this fleld. 
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ABSTRACT 

The design of control laws to damp flexible structure modes requires accurate math 
models of the dynamic system. To obtain the required accuracy of a math model, 
the parameter estimation technique using maximum likelihood estimation is 
employed to improve the accuracy of the model based on flight data. This paper 
presents all phases of this methodology: pre-flight analysis (i.e., optimal input signal 
design for flight test, sensor location determination, model reduction technique, 
etc.), data collection and preprocessing, and post-flight analysis (i.e., estimation 
technique and model verification). The results of this study indicate that the 
parameter estimation technique (i.e., maximum likelihood estimation) is an 
effective and powerful technique in modifing high-order aeroelastic aircraft models. 
However, the accuracy of the results depends upon the fidelity of the theoretical 
model with regards to the correct number of dominant modes for the desired 
frequency bandwith in the model (i.e., model order). If the number of modes in the 
model are not representative, then an identification problem can occure in the 
parameter estimation technique. Nevertheless, this problem can be overcome using 
the system identification technique. 
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INTRODUCTION 


Having an accurate mathematical representation is fundamental to any aircraft 
control system design. In general, aircraft models are developed from a theoretical 
basis and modified by analyzing the experimental data (i.e., wind-tunnel data for 
aerodynamic models or ground shake test data for structural models). Although 
present techniques provide very good dynamic models for the design stages of an 
aircraft, often these models do not match the actual dynamic flight response. This 
problem has generated a need for advanced system identifiaction and parameter 
estimation techniques in upgrading dynamic models of an aircraft based on flight 
test data. This modeling problem is more apparent with high-order aeroelastic 
models with which our experince with modeling techniques is limited. 

Low-frequency structural modes are easily excited for a jet transport with a long 
fuselage. This excitation causes a lateral ride discomfort in certain flight conditions. 
In order to design a yaw damper to dampen Dutch roll response and suppress the 
undesirable low-frequency structure modes by means of active control, an accurate 
aeroelastic model of the aircraft must be available. In this study, parameter 
estimation technique is applied to upgrade the high-order aeroelstic math model of 
a jet transport. The following is a summary of the parameter estimation technique 
using maximum likelihood estimation. 


Maximum Likelihhood Estimation 

Suppose the actual system is described by (Reference 1): 

x(t) = A x(t) + B u(t) + S s (t) + F n (t) 
z(t. ) = C xa.j + Dua.J+Hsfr.J + Gmtt.) 


where 



x(t) 

state vector 

u(t) 

control vector 

z(tj) 

measurement vector 

s (t) 

bias vector 

n(t) 

process noise 

m(t,) 

measurement noise 

ti 

time sample 

A ,B ,C ,D,S ,H,F,G 

system matrices with unknown parameters 


n(t)and m(t) are zero mean Gaussian and independent noise 

Assume k is the vector of unknowns that contains elements of the system matrices 
A, B, C, D, S, H, F and G. The objective is to maximize the probability distribution of 
unknowns (i.e., k) when the measurements z are available. Therefore, maximizing 
P(k/z), where P is the probability distribution function of k given z. 

By Bayes' rule: 


P(k/z) 

P(z) = P(k,z) = P (z / k) P(k) 

(2) 

P(k/z) 

s s 

II 

(3) 


Since in these equations z is given, so P(z) becomes a constant. Assume there is no a 
priori preference for k, so P(k) becomes a constant. Therefore, P(z/k) differs from 
P(k/z) only by a constant. In other words equation (3) becomes: 

P(k/z) = P(z/k) -constant (4) 


Equation (4) indicates that P(z/k) may be maximized instead of P(k/ z). Therefore, 
using Gaussian assumption, the likelihood ratio may be written as: 


I N f N i 

P(z / k) = [(2 b) l IgG'I ] exp|-J £ [*„(!, j-zft,)]* (GGT 1 [z k (t,)-z(t,)]j 


where 


z k (ti) predicted estimate at time ti 

GG* measurement noise covariance matrix 

L number of measurements 


( 5 ) 
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If the logarithm of equation (5) is taken, the consatnt terms are eliminated by the 
maximization, and the equation is multiplied by -1 to do minimization rather than 
maximization, then equation (6) will be obtained as: 

I <k) - 7 f {[z k ft, > - *<',>]* <GGT 1 [z ft > -z ft jj} 

i = 1 (6) 


where J(k) is the cost function to be minimized. Two steps are taken to obtain zjc(tj). 
Prediction step: 

x k^i + i) = ^ x k^i^ + ^ u ^i + l/2^ 

Z k^ t i + l^ = ^' X k^i + 1^ + ^ U ^i + 1^ 

where 

-t 

O = e A * and 'P = J e A * ds 

J o 

and the correction step: 

X k^i + l) = X k ^i + i ) + ^ [ z “ Z k^i + l^] 

K in equation (8) is the Kalman filter gain matrix given by: 

K = P C * (GG* )~ 1 


where P is the solution to the discrete time Riccati equation: 


AP+PA* PC^GG^CP + FF* = 0 


( 10 ) 


After obtaining the cost function J(k), the Newton-Raphson algorithm is used 
iteratively to minimize the cost function by revising the unknowns parameters. 


k 1 ., = k 1 -{V 2 k J(k 1 )}‘'{V k I(k,)} 


( 11 ) 
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This algorithm requires an intial estimate for the vector of unknowns (ko). A priori 
es tim ate is available for each unknown parameter through the analytical model. 

The MMLE software tool developed by NASA Dryden is a parameter estimation 
program supporting this estimation technique. This software has been modified by 
Boeing to accept and handle higher order models. A comprehensive description of 
this software tool is described in Reference 1. 


PRE -FLIGHT ANALYSIS 


Math Model 

A sixtieth order linear aeroelastic math model for a flight condition of Mach .6 
speed, 15000 foot altitude, and no turbulance, and cruise configuration of a jet 
transport was provided in the form of: 


where 


Mq+C 

q + K q = u 

M 

mass matrix 

C 

damping matrix 

K 

stiffness matrix 

q 

generalized coordinate 

u 

control inputs 


( 12 ) 


The model is defined in the inertial axis system, and the dynamics (q), consist of 
rigid body and flexible modes. The model is tuned using data from ground shake 
testing. The system of equations (12) was transformed into state-space form using 
the following transformation: 
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therefore the system equation (12) becomes: 


x = A m x m +B m u 
y = C m x m + D m u 


where 



- M -1 C 


This transformation always exists because the mass matrix is positive definite. 
Although this is a well-posed theoretical problem, it is not trivial. The flexible 
model is usually on the order of one hundred states, thus causing numerical 
inaccuracies in the inversion of the mass matrix. In our analysis the software 
package MPAC was used to perform the transformation. (MPAC is a numerically 
robust modem control and analysis software tool developed by the Boeing 
Company.) 

For the identification process, the system equation (13) was transformed into the 
conjugate modal form using the following transformation: 

m = T -1 x m 


Equation (13) becomes: 

m = Am + B u 
y =C m + D u 

where 
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A = dia (X t ) 
B = T -1 B 
C =CT 


X, = i* eigenvalue 
controllability matrix 
observability matrix 


The advantage of using the modalized form given by equation (14) is that all the 
modes through A matrix, along with the controllability and observability matrices 
are readily available for an analyst to quickly locate uncontrollable and unobservable 
modes. In addition, the modes in the A matrix are decoupled and may be 
partitioned into rigid model and elaastic model. 

The order of the model was reduced to nineteen by deleting the modes above 6 Hz. 
Since this model will eventually be used for ride quality study and modal 
supperasion design, only those modes less than 6 Hz were retained. 

The reduced order, modal model (19th order) is represented by: 

m r = A r m r + B r u 

y =C r m r + D r u (15) 

This model contains one state for heading, one for the spiral mode, two for the 
Dutch roll mode, one for roll mode, eight for low-damped elastic modes, and six for 
high-damped elastic modes. 

To support this study, a special set of sensors were installed on the aircraft to 
measure the dynamic response of the jet transport. The locations of these sensors 
were based on the mode shapes of the aircraft determined by the math model and 
physical constraints (Table I). (A complete discussion on sensor selection and 
location placement on the aircraft is omitted herein for proprietary reasons.) 
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TABLE I: Sensor Type and Locations for High-order 
Aeroelastic Modeling 


SENSOR TYPE 

SENSOR LOCATION 

Position Transducer 

On all control surfaces 

Yaw Rate Gyro 

Pilot seat, IRU (a station between CG and 
cockpit below the cabin floor), CG station 

Lateral Accelerometer 

1 Pilot seat, 1 Cockpit ceiling, 8 on the 
passanger cabin floor from the cockpit to the 
aft galley, 1 on the aft galley ceiling, 3 on 

vertical tail (tip and mid section, front 
and rear spar), three on each nacelle, 

1 IRU station 

Vertical Accelerometer 

1 on the pilot seat, 1 IRU, 1 aft galley, 

8 on each wing, 3 on each horizontal tail, 

2 on each nacelle 

Roll rate, Yaw rate, 
Bank angle, Heading, 

IRU and CG stations 


The sensors selected for the analysis were: body roll angle (O), heading angle OP), roll 
rate (p) and yaw rate (r) at the IRU; body yaw rate at the pilot seat; 9 lateral 
accelerometers along the fuselage; 2 lateral accelerometers on the nacelle number 2; 
and 3 lateral accelerometers on the vertical tail. 


Input Signal Desig n 

The flight test input-signal design analysis for high-order aeroelastic modeling was 
performed using the reduced order analytical model (equation 15). Although a 
number of "optimum" input signals have been proposed for flight testing in 
conjunction with parameter estimation, none have been found to be appropriate for 
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high-order aeroelastic modeling. Essentially, all the analytical techniques proposed 
in designing the optimum input signals are based on the analytical model. This 
model is the subject of improvement by the identification and estimation 
techniques. Hence, no "optimum" input signal exists. 

A number of different input signals were evaluated for this study. After a 
comprehensive simulation study, it was determined that a frequency sweep of a 
Hnoar «:inp-wavp with adequate energy to excite all the modes (rigid and elastic) 
yeilds the best results. In addition, the linear sine-wave frequency sweep optimizes 
the most commonly used criterion for input signal design: 

91 = - log (det M) (16) 

where M is the Fisher information matrix (or sensitivity matrix) defined by: 

M = V*J(k) (17) 

I is (he cost function defined in equation (6). The criterion 9i defined in equation 

(16) is related to the volume of highest probability density region for the parameters 
k. An interesting property of the determinant criterion is that it is independent of 
scaling parameters (Refemce 2). 

Fifteen tests were designed for the same flight condition. Five frequency sweeps 
were designed for each control surface. Each test was repeated for rudder, aileron, 
and both surfaces in phase. The first frequency sweep covered 0 to 6 Hz to excite all 
the modes in one test. The other four tests were then designed to excite specifically 
high-damped modes by sweeping from .25 Hz below to .25 Hz above the frequency of 
the mode. 

The amplitude of the input signals were designed to be constant for practical 
purposes (i.e., rate limits). The designed input signals were tested in the lab to 
confirm that the signals did not saturate the servos and actuators of the control 
surfaces. However, the output of the actuators during flight test generated signals 
with decaying amplitutes. These decaying amplitudes reduced the energy level 
initially designed for the test. Figures 1 and 2 show the actual control surface 
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Figure 1. Actual Control Surface Deflections for Flight Condition 21 





Figure 2. Actual Control Surface Deflections for Flight Condition 41 



deflections for rudder sweep alone, and for rudder and aileron surfaces 
simultaneously in phase. 


Sampling Frequency 

To record the data in flight test, a simulation study was conducted to determine the 
required sampling frequency. The analytical model (i.e., system equations 15) was 
assumed to be the true model, and simulated using the designed input sig nal . A 
considerable amount of noise was added to the simulation data, and then that data 
was treated as pseudo-flight data. The acutal model was used for parameter 
estimation to determine the required sampling frequency. Sampling frequencies of 
20, 25, 50, 100, 200 Hz were considered for this study. One mode or group of modes 
at a time were selected for the estimation process of each sampling frequency. The 
results indicated that 100 Hz is the best sampling frequency for this study. Figure 3 
shows the typical results for identified parameters when different sampling 
frequencies were used. 



Figure 3. Typical Results from Estimation with Different 
Sampling Frequency 
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FLIGHTIEST 


The flight test was performed using designed linear sine-wave frequency sweeps for 
rudder and aileron. The test conditions were conducted at a speed of Mach .6, an 
altitude of 15,000 feet, and minimal turbulence. A preprogrammed frequency 
function generator was used to apply the linear sinusoidal frequency sweeps (0-6 Hz) 
to the aileron and rudder (through the autopilot servo). 

The flight test data were recorded with 100 sample per second, and then filtered 
using a Graham low-pass filter with the cutoff frequency of 10 Hz and rolloff 
frequency of 15 Hz. Prior to estimation analysis, the data were cleaned up by 
removing all the sensor biases and data dropouts. 


POST FLIGHT ANALYSIS 

The analytical model (system equations 15) was simulated using actual control 
surface defelection during flight as input signals. The comparison of flight data with 
the response of the analytical model for flight condition 41, where both rudder and 
aileron frequency sweeps are used, is presented in the Figures 4-11. 

The maximum likelihood estimation software tool (MMLE) developed by NASA 
Dryden was used to minimize the residuals between flight data and response of the 
analytical model in Figures 4-11. At the time of analysis, MMLE was hosted on the 
Cyber mainfram. Due to Cyber having a memory limit, the capability of using 
process noise was not available for analysis. Hence the results obtained herein, are 
preliminary results which do not include the effect of process noise. The final 
results of this study will be reported at the 1989 AIAA Guidance, Navigation and 
Control conference. 

The high-order model was partitioned into two sections: rigid model and elastic 
model. For rigid model identification, 15 seconds of data were used. First the rigid 
portion of the control and measurement matrices were upgraded. Then, the A 
matrix was upgraded. Finally, all the parameters in the rigid section of the A /B and 
C matrices were simultaneously estimated. 
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Flight condition 41, Rudder and Aileron Input 



Figure 4. Time Response of Body Roll at IRU Comparing Flight Data with Math Model 



Flight condition 41, Ruddor and AJlwon Input 



Figure 5. Time Response of Body Yaw at IRU Comparing Flight Data with Math Model 









Right condition 41, Rudder and Aileron Input 



Frequency, HZ 

Figure 7. Time Response of Yaw Rate at IRU Comparing Flight Data with Math Model 
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Right condition 41, Rudder and Aileron input 



Frequency, HZ 

Figure 10. Time Response of Lateral Acceleration at AFT Body Comparing Flight Data 
with Math Model 




Flight condition 41, Rudder and Aileron input 



Figure 11. Time Response of Lateral Acceleration at Engine 2 Inlet Comparing Flight Data 
with Math Model 





Two different approches were taken for the elastic model identification. In the first 
approach, die 19 th order model was used for the analysis with all the elements of B 
and C being estimated. All 70 seconds of data were used for this estimation 
approach, . In this process, those parameters in the B andC that did not contribute 
to the residuals were identified and kept constant for the remainder of the analysis. 
Then, the elements of A were added to the estimation process while keeping some 
of the elements of B andC constant. The results of this estimation approach are 
show in Figures 12-19. 

The second approach was to add one elastic mode at a time to the rigid model. For 
this approach, the first elastic mode was added with 28 seconds of data used for the 
analysis. The corresponding parameters in the B and C matrices were estimated 
every time a mode was added to the model. The result of this approach was not 
satisfactory because several times the algorithem diverged and the residuals were 
big. 

Figures 20-26 show the PSD plots obtained from the analytical model. Figures 27-34 
show the PSD plots obtained from the estimated model. The PSD plots obtained 
from the estimated model, clearly show that the estimation analysis improved the 
accuracy of the model in terms of its modal representation. However, the estimated 
parameters in the B and C matrices are biased. Since an accurate representation of 
the transfer functions was desired for this study rather than true values of the B 
and C matrices, the biased estimates in theB and C matrices did not create any 
problem. 

Figures 16, 17 and 19 indicate that another mode is present in the flight data which is 
not modeled in the analytical or estimated model. This problem can not be solved 
via parameter estimation technique which assumes the structure of the model (i.e., 
the order of the model) is correct. Hence, it is suggested that the system 
identification technique developed by V. klein and J. Batterson of NASA LaRC be 
used to overcome this problem. 
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Flight condition 41, Rudder and Aileron Input 
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Frequency, HZ 

Figure 12. Time Response of Body Roll at IRU Comparing Flight Data with MLE Model 





Right condition 41, Rudder and Aileron Input 



Figure 13. Time Response of Body Yaw at IRU Comparing Flight Data with MLE Model 




Flight condition 41, Rudder and Aileron Input 



Figure 15. Time Response of Yaw Rate at IRU Comparing Flight Data with MLE Model 
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Figure 20. Power Spectral Density Response of Body Roll at IRU Comparing Flight Data 
with Math Model 
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Figure 21. Power Spectral Density Response of Body Yaw at IRU Comparing Flight Data 
with Math Model 
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Figure 22. Power Spectral Density Response of Roll Rate at IRU Comparing Flight Data 
with Math Model 
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Frequency, HZ 

Figure 23. Power Spectral Density Response of Yaw Rate at IRU Comparing Flight Data 
with Math Model 
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Figure 24. Power Spectral Density Response of Lateral Acceleration at Pilot Seat Comparing 
Flight Data with Math Model 
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9 *o. Power Spectral Density Response of Lateral Acceleration at IRU Comparing 

Flight Data with Math Model 
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Figure 26. Power Spectral Density Response of Lateral Acceleration at AFT Body 
Comparing Flight Data with Math Model 
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Figure 28. Power Spectral Density Response of Body Yaw at IRU Comparing Flight Data 
with MLE Model 
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Figure 29 . Power Spectral Density Response of Roll Rate at IRU Comparing Flight Data 
with MLE Model 
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Figure 31. Power Spectral Density Response of Lateral Acceleration at Pilot Seat 
Comparing Flight Data with NILE Model 
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Figure 32. Power Spectral Density Response of Lateral Acceleration at IRU 
Comparing Flight Data with MLE Model 
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Figure 33. Power Spectral Density Response of Lateral Acceleration at AFT Body 
Comparing Flight Data with MLE Model 
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Frequency, HZ 

Figure 34. Power Spectral Density of Lateral Acceleration at Engine 2 Inlet 
Comparing Flight Data with WILE Model 
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